"""
Comprehensive Test Suite for Thermal-Fluid Physics

Tests for all thermal and multi-physics capabilities in NSEMSolver.jl:
1. Material Properties Framework
2. Boussinesq Natural Convection  
3. Heat Transfer Solvers
4. Multi-Physics Coupling
5. Thermal Boundary Conditions
6. Turbulent Heat Transfer
7. Thermal Visualization
8. Validation Against Analytical Solutions

Ensures research-grade accuracy and reliability.
"""

using Test
using NSEMSolver
using LinearAlgebra
using Statistics

# Test tolerances
const RTOL_LOOSE = 1e-2   # For convergence tests
const RTOL_TIGHT = 1e-4   # For property tests
const RTOL_STRICT = 1e-6  # For analytical comparisons

@testset "Thermal-Fluid Physics Tests" begin

@testset "Material Properties Framework" begin
    @testset "Constant Properties" begin
        # Test constant property model
        props = ConstantProperties(
            density = 1.177,
            viscosity = 1.85e-5, 
            thermal_conductivity = 0.026,
            specific_heat = 1007.0
        )
        
        @test compute_density(300.0, 101325.0, props) ≈ 1.177 rtol=RTOL_STRICT
        @test compute_viscosity(350.0, 101325.0, props) ≈ 1.85e-5 rtol=RTOL_STRICT
        @test compute_thermal_conductivity(250.0, 101325.0, props) ≈ 0.026 rtol=RTOL_STRICT
        @test compute_specific_heat(400.0, 101325.0, props) ≈ 1007.0 rtol=RTOL_STRICT
        
        # Test derived properties
        α = compute_thermal_diffusivity(300.0, 101325.0, props)
        @test α ≈ 0.026 / (1.177 * 1007.0) rtol=RTOL_TIGHT
        
        ν = compute_kinematic_viscosity(300.0, 101325.0, props)  
        @test ν ≈ 1.85e-5 / 1.177 rtol=RTOL_TIGHT
        
        Pr = compute_prandtl_number(300.0, 101325.0, props)
        @test Pr ≈ ν / α rtol=RTOL_TIGHT
    end
    
    @testset "Temperature-Dependent Properties" begin
        # Test air properties with temperature dependence
        air_props = create_air_properties(model=:temperature_dependent)
        
        # Properties should vary with temperature
        T1, T2 = 250.0, 350.0
        P = 101325.0
        
        μ1 = compute_viscosity(T1, P, air_props)
        μ2 = compute_viscosity(T2, P, air_props)
        @test μ2 > μ1  # Viscosity increases with temperature for gases
        
        k1 = compute_thermal_conductivity(T1, P, air_props)
        k2 = compute_thermal_conductivity(T2, P, air_props)
        @test k2 > k1  # Thermal conductivity increases with temperature
        
        ρ1 = compute_density(T1, P, air_props)
        ρ2 = compute_density(T2, P, air_props)
        @test ρ2 < ρ1  # Density decreases with temperature (ideal gas)
        
        # Test realistic ranges
        @test 200.0 < μ1 * 1e7 < 300.0  # Viscosity in reasonable range
        @test 0.02 < k1 < 0.04          # Thermal conductivity range
        @test 1.0 < ρ1 < 2.0            # Density range
    end
    
    @testset "Sutherland Viscosity Model" begin
        sutherland = SutherlandModel(
            mu0 = 1.716e-5,
            T0 = 273.15,
            S = 110.4
        )
        
        # Test Sutherland's law at reference conditions
        μ_ref = compute_viscosity(273.15, 101325.0, sutherland)
        @test μ_ref ≈ 1.716e-5 rtol=RTOL_TIGHT
        
        # Test temperature dependence  
        μ_300 = compute_viscosity(300.0, 101325.0, sutherland)
        μ_400 = compute_viscosity(400.0, 101325.0, sutherland)
        @test μ_400 > μ_300 > μ_ref
        
        # Test against known values (approximate)
        μ_300_expected = 1.85e-5  # Approximate viscosity at 300K
        @test μ_300 ≈ μ_300_expected rtol=0.1
    end
    
    @testset "Fluid Property Validation" begin
        # Test air properties creation
        air = create_air_properties()
        water = create_water_properties()
        
        # Basic validation tests
        @test validate_properties(air, 300.0, 101325.0)
        @test validate_properties(water, 300.0, 101325.0)
        
        # Test property summary (should not error)
        @test_nowarn print_property_summary(air, 300.0, 101325.0)
    end
end

@testset "Boussinesq Natural Convection" begin
    @testset "Boussinesq Options Construction" begin
        opts = BoussinesqOptions(
            reference_temperature = 300.0,
            thermal_expansion = 1/300.0,
            gravitational_acceleration = [0.0, -9.81],
            density_reference = 1.177,
            prandtl_number = 0.71
        )
        
        @test opts.reference_temperature == 300.0
        @test opts.thermal_expansion ≈ 1/300.0
        @test opts.gravitational_acceleration == [0.0, -9.81]
        @test opts.prandtl_number == 0.71
    end
    
    @testset "Rayleigh Number Computation" begin
        bouss_opts = BoussinesqOptions(
            reference_temperature = 300.0,
            thermal_expansion = 1/300.0,
            gravitational_acceleration = [0.0, -9.81]
        )
        
        ΔT = 10.0
        L = 1.0
        Ra = compute_rayleigh_number(bouss_opts, ΔT, L)
        
        # Should be positive and reasonable
        @test Ra > 0
        @test 1e3 < Ra < 1e8  # Typical range
        
        # Test scaling: Ra ~ ΔT
        Ra2 = compute_rayleigh_number(bouss_opts, 2*ΔT, L)
        @test Ra2 ≈ 2*Ra rtol=RTOL_TIGHT
        
        # Test scaling: Ra ~ L³
        Ra3 = compute_rayleigh_number(bouss_opts, ΔT, 2*L)
        @test Ra3 ≈ 8*Ra rtol=RTOL_TIGHT
    end
    
    @testset "Buoyancy Force Computation" begin
        nx, ny = 5, 5
        T = 300.0 .+ 10.0 * rand(nx, ny)  # Temperature field
        buoyancy_force = zeros(nx, ny, 2)
        
        ρ_ref = 1.177
        β = 1/300.0  
        g_vec = [0.0, -9.81]
        T_ref = 300.0
        
        compute_buoyancy_force!(buoyancy_force, T, ρ_ref, β, g_vec, T_ref)
        
        # Test dimensions
        @test size(buoyancy_force) == (nx, ny, 2)
        
        # Test physics: warmer fluid should have upward buoyancy
        for i in 1:nx, j in 1:ny
            ΔT = T[i,j] - T_ref
            expected_fy = -ρ_ref * β * ΔT * g_vec[2]  # Negative g, so hot fluid goes up
            @test buoyancy_force[i,j,2] ≈ expected_fy rtol=RTOL_TIGHT
        end
        
        # Horizontal component should be zero for vertical gravity
        @test all(abs.(buoyancy_force[:,:,1]) .< 1e-10)
    end
    
    @testset "Simple Natural Convection Solve" begin
        # Test minimal natural convection case
        ns_opts = NSOptions(N=3, nu=1e-4, cfl=0.5, tfinal=1.0, dim=2, 
                           n_block=2, max_steps=50, tol=1e-2, verbose=false)
        
        bouss_opts = BoussinesqOptions(
            reference_temperature = 300.0,
            thermal_expansion = 1/300.0, 
            gravitational_acceleration = [0.0, -9.81],
            prandtl_number = 0.71
        )
        
        # This should not error and produce a result
        @test_nowarn result = solve_boussinesq_2d(ns_opts, bouss_opts)
        
        result = solve_boussinesq_2d(ns_opts, bouss_opts)
        
        # Basic result validation
        @test !isnothing(result.u)
        @test !isnothing(result.v) 
        @test !isnothing(result.T)
        @test size(result.u) == size(result.v) == size(result.T)
        @test length(result.x) == size(result.u, 1)
        @test length(result.y) == size(result.u, 2)
    end
end

@testset "Heat Transfer Solvers" begin
    @testset "Heat Transfer Options" begin
        heat_opts = HeatTransferOptions(
            thermal_conductivity = 0.026,
            thermal_diffusivity = 2.2e-5,
            heat_capacity = 1007.0,
            density = 1.177,
            boundary_conditions = Dict{Symbol,Any}(
                :bottom => Dict(:type => :dirichlet, :value => 310.0),
                :top => Dict(:type => :dirichlet, :value => 290.0),
                :left => Dict(:type => :neumann, :value => 0.0),
                :right => Dict(:type => :neumann, :value => 0.0)
            )
        )
        
        @test heat_opts.thermal_conductivity == 0.026
        @test heat_opts.thermal_diffusivity ≈ 2.2e-5
        @test haskey(heat_opts.boundary_conditions, :bottom)
        @test heat_opts.boundary_conditions[:bottom][:type] == :dirichlet
    end
    
    @testset "Heat Equation Solution" begin
        # Test basic heat equation solver
        nx, ny = 6, 6
        x = LinRange(-1, 1, nx)
        y = LinRange(-1, 1, ny)
        
        T = zeros(nx, ny)
        T_old = 300.0 .+ 10.0 * sin.(π * x') .* sin.(π * y)  # Smooth initial condition
        
        u = zeros(nx, ny)  # No convection
        v = zeros(nx, ny)
        
        dt = 1e-4
        α = 2.2e-5
        
        # Apply heat equation step
        solve_heat_equation_2d!(T, T_old, u, v, dt, α, x, y)
        
        # Temperature should change due to diffusion
        @test T != T_old
        
        # Total heat should be approximately conserved (no sources/sinks)
        Q_old = sum(T_old)
        Q_new = sum(T)
        @test Q_new ≈ Q_old rtol=RTOL_LOOSE  # Some boundary effects expected
        
        # Should be stable (no blow-up)
        @test all(isfinite.(T))
        @test maximum(abs.(T)) < 2 * maximum(abs.(T_old))
    end
    
    @testset "Nusselt Number Computation" begin
        # Create simple temperature field with known gradient
        nx, ny = 10, 10
        x = LinRange(-1, 1, nx)
        y = LinRange(-1, 1, ny)
        
        # Linear temperature profile: T = T0 + (y+1) * ΔT/2
        T0, ΔT = 290.0, 20.0
        T = T0 .+ (y .+ 1)' .* ΔT/2
        
        L = 2.0  # Domain height
        k = 0.026  # Thermal conductivity
        
        Nu = compute_nusselt_number_2d(T, :bottom, L, k)
        
        # For linear profile, Nu should be close to 1 (pure conduction)
        @test Nu ≈ 1.0 rtol=RTOL_LOOSE
        
        # Test different wall
        Nu_top = compute_nusselt_number_2d(T, :top, L, k)
        @test Nu_top ≈ 1.0 rtol=RTOL_LOOSE
    end
end

@testset "Multi-Physics Coupling" begin
    @testset "Multi-Physics Options" begin
        mp_opts = MultiPhysicsOptions(
            include_heat_transfer = true,
            include_boussinesq = true,
            coupling_strategy = :segregated,
            max_coupling_iterations = 10,
            coupling_tolerance = 1e-5
        )
        
        @test mp_opts.include_heat_transfer == true
        @test mp_opts.include_boussinesq == true
        @test mp_opts.coupling_strategy == :segregated
        @test mp_opts.max_coupling_iterations == 10
    end
    
    @testset "Coupled System Time Step" begin
        ns_opts = NSOptions(nu=1e-4, cfl=0.3)
        mp_opts = MultiPhysicsOptions(include_heat_transfer=true, include_turbulence=false)
        
        x = LinRange(-1, 1, 20)
        y = LinRange(-1, 1, 20)
        
        dt = compute_multiphysics_timestep_2d(ns_opts, mp_opts, x, y)
        
        # Should be positive and reasonable
        @test dt > 0
        @test dt < 1.0  # Should be small for stability
        
        # Should be smaller with more physics included
        mp_opts_simple = MultiPhysicsOptions(include_heat_transfer=false)
        dt_simple = compute_multiphysics_timestep_2d(ns_opts, mp_opts_simple, x, y)
        @test dt ≤ dt_simple
    end
    
    @testset "Coupled Residual Computation" begin
        nx, ny = 5, 5
        u = randn(nx, ny) * 0.1
        v = randn(nx, ny) * 0.1  
        T = 300.0 .+ randn(nx, ny) * 5.0
        
        residual = compute_coupled_residual_2d(u, v, T)
        
        @test residual ≥ 0
        @test isfinite(residual)
        
        # Zero fields should give zero residual
        u_zero = zeros(nx, ny)
        v_zero = zeros(nx, ny)
        T_zero = zeros(nx, ny)
        
        residual_zero = compute_coupled_residual_2d(u_zero, v_zero, T_zero)
        @test residual_zero ≈ 0.0 rtol=RTOL_STRICT
    end
    
    @testset "Simple Coupled Flow-Heat Solve" begin
        # Minimal coupled simulation test
        ns_opts = NSOptions(N=3, nu=1e-4, cfl=0.8, tfinal=0.5, dim=2,
                           n_block=2, max_steps=20, tol=1e-1, verbose=false)
        
        bouss_opts = BoussinesqOptions(
            reference_temperature = 300.0,
            thermal_expansion = 1/300.0,
            gravitational_acceleration = [0.0, -9.81],
            prandtl_number = 0.71
        )
        
        heat_opts = HeatTransferOptions(
            thermal_conductivity = 0.026,
            thermal_diffusivity = 2.2e-5,
            boundary_conditions = Dict{Symbol,Any}(
                :bottom => Dict(:type => :dirichlet, :value => 310.0),
                :top => Dict(:type => :dirichlet, :value => 290.0),
                :left => Dict(:type => :neumann, :value => 0.0),
                :right => Dict(:type => :neumann, :value => 0.0)
            )
        )
        
        mp_opts = MultiPhysicsOptions(
            coupling_strategy = :segregated,
            max_coupling_iterations = 5,
            coupling_tolerance = 1e-3
        )
        
        @test_nowarn result = solve_coupled_flow_heat_2d(ns_opts, bouss_opts, heat_opts, mp_opts)
        
        result = solve_coupled_flow_heat_2d(ns_opts, bouss_opts, heat_opts, mp_opts)
        
        # Validate result structure
        @test result isa MultiPhysicsResult
        @test !isnothing(result.u) && !isnothing(result.v)
        @test !isnothing(result.T)
        @test !isnothing(result.physics_diagnostics)
        @test haskey(result.physics_diagnostics, "average_temperature")
        @test result.solve_time ≥ 0
    end
end

@testset "Thermal Boundary Conditions" begin
    @testset "Dirichlet Temperature BC" begin
        dirichlet_bc = DirichletTemperature(temperature=350.0)
        @test dirichlet_bc.temperature == 350.0
        @test dirichlet_bc.time_dependent == false
        
        # Time-dependent case
        temp_func = (x, y, t) -> 300.0 + 10.0 * sin(t)
        dirichlet_td = DirichletTemperature(temperature=temp_func, time_dependent=true)
        @test dirichlet_td.time_dependent == true
        @test dirichlet_td.temperature(0, 0, π/2) ≈ 310.0 rtol=RTOL_TIGHT
    end
    
    @testset "Neumann Heat Flux BC" begin
        neumann_bc = NeumannHeatFlux(heat_flux=1000.0)  # W/m²
        @test neumann_bc.heat_flux == 1000.0
        @test neumann_bc.time_dependent == false
        
        # Variable heat flux
        flux_func = (x, y, t) -> 1000.0 * (1 + 0.1 * cos(2π * t))
        neumann_var = NeumannHeatFlux(heat_flux=flux_func, time_dependent=true)
        @test neumann_var.time_dependent == true
        @test neumann_var.heat_flux(0, 0, 0) ≈ 1100.0 rtol=RTOL_TIGHT
    end
    
    @testset "Convective Boundary BC" begin
        conv_bc = ConvectiveBoundary(
            heat_transfer_coefficient = 25.0,  # W/m²·K
            ambient_temperature = 280.0        # K
        )
        
        @test conv_bc.heat_transfer_coefficient == 25.0
        @test conv_bc.ambient_temperature == 280.0
        @test conv_bc.time_dependent == false
    end
    
    @testset "Thermal Boundary System" begin
        boundary_system = create_thermal_boundary_system(thermal_conductivity=0.026)
        @test boundary_system.thermal_conductivity == 0.026
        
        # Add boundary conditions
        hot_wall = DirichletTemperature(temperature=350.0)
        cold_wall = DirichletTemperature(temperature=250.0)
        
        add_boundary_condition!(boundary_system, :left, hot_wall)
        add_boundary_condition!(boundary_system, :right, cold_wall)
        
        @test haskey(boundary_system.boundary_conditions, :left)
        @test haskey(boundary_system.boundary_conditions, :right)
        
        # Time update
        update_time!(boundary_system, 1.0, 0.01)
        @test boundary_system.current_time == 1.0
        @test boundary_system.time_step == 0.01
    end
    
    @testset "Boundary Condition Application" begin
        # Create simple test case
        nx, ny = 6, 6
        T = zeros(nx, ny) .+ 300.0  # Initial temperature
        x = LinRange(-1, 1, nx)
        y = LinRange(-1, 1, ny)
        
        # Create boundary system
        boundary_system = create_thermal_boundary_system()
        
        # Add Dirichlet conditions
        hot_bc = DirichletTemperature(temperature=350.0)
        cold_bc = DirichletTemperature(temperature=250.0)
        
        add_boundary_condition!(boundary_system, :left, hot_bc)
        add_boundary_condition!(boundary_system, :right, cold_bc)
        
        # Apply boundary conditions
        apply_thermal_boundary_conditions!(T, boundary_system, x, y)
        
        # Check that boundaries were applied
        @test all(T[1, :] .≈ 350.0)    # Left wall hot
        @test all(T[end, :] .≈ 250.0)  # Right wall cold
        
        # Interior should be unchanged initially
        @test T[2, 2] ≈ 300.0
    end
end

@testset "Turbulence Models" begin
    @testset "Smagorinsky Model" begin
        smag_model = SmagorinskyModel(C_s=0.17, filter_width=1.0)
        @test smag_model.C_s == 0.17
        @test smag_model.filter_width == 1.0
        @test smag_model.van_driest_damping == false
    end
    
    @testset "Dynamic Smagorinsky Model" begin
        dyn_smag = DynamicSmagorinskyModel(
            test_filter_ratio = 2.0,
            clip_coefficient = true,
            minimum_coefficient = 0.0
        )
        
        @test dyn_smag.test_filter_ratio == 2.0
        @test dyn_smag.clip_coefficient == true
        @test dyn_smag.minimum_coefficient == 0.0
    end
    
    @testset "WALE Model" begin
        wale_model = WALEModel(C_w=0.325, filter_width=1.0)
        @test wale_model.C_w == 0.325
        @test wale_model.filter_width == 1.0
    end
    
    @testset "Strain Rate Tensor Computation" begin
        nx, ny = 5, 5
        x = LinRange(-1, 1, nx)
        y = LinRange(-1, 1, ny)
        
        # Simple shear flow: u = y, v = 0
        u = zeros(nx, ny)
        v = zeros(nx, ny)
        for (i, xi) in enumerate(x), (j, yj) in enumerate(y)
            u[i, j] = yj  # Linear shear
        end
        
        strain_rate = zeros(nx, ny)
        compute_strain_rate_tensor_2d!(strain_rate, u, v, x, y)
        
        # For simple shear u = y, strain rate should be constant (≈ 1)
        interior_strain = strain_rate[2:end-1, 2:end-1]
        @test all(interior_strain .> 0)  # Positive strain rate
        @test maximum(interior_strain) < 3.0  # Reasonable magnitude
        
        # Check dimensions
        @test size(strain_rate) == (nx, ny)
    end
    
    @testset "Subgrid Viscosity Computation" begin
        nx, ny = 5, 5
        x = LinRange(-1, 1, nx)  
        y = LinRange(-1, 1, ny)
        
        u = randn(nx, ny) * 0.1
        v = randn(nx, ny) * 0.1
        strain_rate = abs.(randn(nx, ny)) .+ 0.1  # Positive strain rates
        
        nu_sgs = zeros(nx, ny)
        
        # Test Smagorinsky model
        smag_model = SmagorinskyModel(C_s=0.17)
        compute_subgrid_viscosity_2d!(nu_sgs, u, v, strain_rate, smag_model, x, y)
        
        @test all(nu_sgs .≥ 0)  # Non-negative viscosity
        @test size(nu_sgs) == (nx, ny)
        @test maximum(nu_sgs) > 0  # Should produce some subgrid viscosity
        @test all(isfinite.(nu_sgs))
    end
    
    @testset "Explicit Filtering" begin
        nx, ny = 6, 6
        field = randn(nx, ny)
        
        # Test top-hat filter
        filtered_tophat = apply_explicit_filter_2d(field, :top_hat)
        @test size(filtered_tophat) == size(field)
        @test filtered_tophat != field  # Should change the field
        
        # Filtered field should be smoother (lower variance for random input)
        @test var(filtered_tophat) ≤ var(field)
        
        # Test Gaussian filter
        filtered_gauss = apply_explicit_filter_2d(field, :gaussian)
        @test size(filtered_gauss) == size(field)
        @test filtered_gauss != field
        
        # Test filter preservation of mean (approximately)
        @test abs(mean(filtered_tophat) - mean(field)) < 0.5 * std(field)
    end
end

@testset "Validation Against Analytical Solutions" begin
    @testset "Pure Conduction Validation" begin
        # Test against 1D steady-state conduction: T(x) = T0 + (T1-T0)*x
        # Between two walls at different temperatures
        
        nx, ny = 21, 11  
        x = LinRange(0, 1, nx)  # 0 to 1 domain
        y = LinRange(-0.5, 0.5, ny)
        
        T0, T1 = 300.0, 400.0  # Wall temperatures
        
        # Analytical solution: linear profile
        T_analytical = zeros(nx, ny)
        for (i, xi) in enumerate(x)
            T_analytical[i, :] .= T0 + (T1 - T0) * xi
        end
        
        # Set up heat transfer problem
        heat_opts = HeatTransferOptions(
            thermal_conductivity = 1.0,
            thermal_diffusivity = 1.0,
            boundary_conditions = Dict{Symbol,Any}(
                :left => Dict(:type => :dirichlet, :value => T0),
                :right => Dict(:type => :dirichlet, :value => T1),
                :bottom => Dict(:type => :neumann, :value => 0.0),
                :top => Dict(:type => :neumann, :value => 0.0)
            )
        )
        
        # Solve steady-state (run until convergence)
        T = zeros(nx, ny) .+ (T0 + T1)/2  # Initial guess
        
        # Simple iterative solver for steady-state
        for iter in 1:1000
            T_old = copy(T)
            
            # Apply diffusion step
            dt = 1e-3
            α = 1.0
            u = zeros(nx, ny)  # No convection
            v = zeros(nx, ny)
            
            solve_heat_equation_2d!(T, T_old, u, v, dt, α, x, y)
            
            # Apply boundary conditions
            T[1, :] .= T0   # Left wall
            T[end, :] .= T1 # Right wall
            
            # Check convergence
            if maximum(abs.(T - T_old)) < 1e-6
                break
            end
        end
        
        # Compare with analytical solution
        # Interior points should match linear profile
        T_centerline = T[:, div(ny, 2)]
        T_analytical_centerline = T_analytical[:, div(ny, 2)]
        
        # Allow some tolerance for numerical solution
        @test maximum(abs.(T_centerline - T_analytical_centerline)) < 0.01 * (T1 - T0)
        
        # Test heat flux constancy
        heat_flux = -(T[2:end, div(ny,2)] - T[1:end-1, div(ny,2)]) / (x[2] - x[1])
        @test std(heat_flux) / mean(abs.(heat_flux)) < 0.05  # Heat flux should be constant
    end
    
    @testset "Natural Convection Scaling Laws" begin
        # Test that Nusselt number scales appropriately with Rayleigh number
        # Nu ~ Ra^0.25 for laminar natural convection
        
        Ra_values = [1e4, 2e4, 4e4]  # Test range
        Nu_values = Float64[]
        
        for Ra in Ra_values
            # Set up minimal natural convection problem
            ns_opts = NSOptions(N=3, nu=1e-4, cfl=0.5, tfinal=0.5, dim=2,
                               n_block=2, max_steps=10, tol=1e-1, verbose=false)
            
            # Compute properties to achieve target Ra
            L = 1.0
            ΔT = 10.0
            β = 1/300.0
            g = 9.81
            Pr = 0.71
            
            ν = sqrt(g * β * ΔT * L^3 / (Ra * Pr))
            ns_opts = NSOptions(ns_opts; nu=ν)  # Update viscosity
            
            bouss_opts = BoussinesqOptions(
                reference_temperature = 300.0,
                thermal_expansion = β,
                gravitational_acceleration = [0.0, -g],
                prandtl_number = Pr
            )
            
            heat_opts = HeatTransferOptions(
                thermal_conductivity = 0.026,
                thermal_diffusivity = ν / Pr,
                boundary_conditions = Dict{Symbol,Any}(
                    :bottom => Dict(:type => :dirichlet, :value => 305.0),
                    :top => Dict(:type => :dirichlet, :value => 295.0),
                    :left => Dict(:type => :neumann, :value => 0.0),
                    :right => Dict(:type => :neumann, :value => 0.0)
                )
            )
            
            mp_opts = MultiPhysicsOptions(coupling_tolerance=1e-2)
            
            try
                result = solve_coupled_flow_heat_2d(ns_opts, bouss_opts, heat_opts, mp_opts)
                Nu = get(result.physics_diagnostics, "nusselt_number", NaN)
                if isfinite(Nu) && Nu > 0
                    push!(Nu_values, Nu)
                else
                    push!(Nu_values, NaN)  # Mark as failed
                end
            catch
                push!(Nu_values, NaN)  # Solver failed
            end
        end
        
        # Check that we got some valid results
        valid_results = sum(.!isnan.(Nu_values))
        @test valid_results ≥ 2  # At least 2 successful cases
        
        if valid_results ≥ 2
            # Check basic scaling trend (Nu should increase with Ra)
            valid_indices = .!isnan.(Nu_values)
            Ra_valid = Ra_values[valid_indices]
            Nu_valid = Nu_values[valid_indices]
            
            if length(Nu_valid) ≥ 2
                # Nu should generally increase with Ra
                correlation = cor(log.(Ra_valid), log.(Nu_valid))
                @test correlation > 0.3  # Positive correlation (loose tolerance)
            end
        end
    end
end

@testset "Error Handling and Edge Cases" begin
    @testset "Invalid Parameters" begin
        # Test invalid Rayleigh number computation
        bouss_opts = BoussinesqOptions(thermal_expansion = 0.0)  # Invalid β
        @test_throws Exception compute_rayleigh_number(bouss_opts, 10.0, 1.0)
        
        # Test zero temperature difference  
        valid_bouss = BoussinesqOptions()
        @test compute_rayleigh_number(valid_bouss, 0.0, 1.0) ≈ 0.0
        
        # Test negative length
        @test_throws Exception compute_rayleigh_number(valid_bouss, 10.0, -1.0)
    end
    
    @testset "Boundary Condition Errors" begin
        # Test missing boundary conditions
        incomplete_bc = Dict{Symbol,Any}(
            :bottom => Dict(:type => :dirichlet, :value => 300.0)
            # Missing other walls
        )
        
        heat_opts = HeatTransferOptions(boundary_conditions = incomplete_bc)
        
        # Should handle missing BCs gracefully (use defaults)
        @test_nowarn heat_opts
    end
    
    @testset "Solver Robustness" begin
        # Test with extreme parameters
        ns_opts = NSOptions(N=2, nu=1e-8, max_steps=5, tol=1e-1, verbose=false)
        
        bouss_opts = BoussinesqOptions(
            thermal_expansion = 1e-2,  # Very large expansion
            gravitational_acceleration = [0.0, -100.0]  # Very large gravity
        )
        
        # This might fail, but should not crash
        @test_nowarn try
            solve_boussinesq_2d(ns_opts, bouss_opts)
        catch e
            # Expected to fail with extreme parameters
        end
    end
end

@testset "Performance and Memory" begin
    @testset "Memory Allocation" begin
        nx, ny = 10, 10
        x = LinRange(-1, 1, nx)
        y = LinRange(-1, 1, ny)
        
        # Test that functions don't allocate excessively
        u = randn(nx, ny) * 0.1
        v = randn(nx, ny) * 0.1
        strain_rate = zeros(nx, ny)
        
        # This should not allocate much memory
        @test (@allocated compute_strain_rate_tensor_2d!(strain_rate, u, v, x, y)) < 1000
        
        # Buoyancy force computation
        T = zeros(nx, ny) .+ 300.0
        buoyancy_force = zeros(nx, ny, 2)
        
        @test (@allocated compute_buoyancy_force!(buoyancy_force, T, 1.0, 1/300.0, [0.0, -9.81], 300.0)) < 1000
    end
    
    @testset "Scaling Behavior" begin
        # Test that solvers scale reasonably with problem size
        function time_solve(n_block)
            ns_opts = NSOptions(N=3, nu=1e-4, tfinal=0.1, dim=2, n_block=n_block,
                               max_steps=5, tol=1e-1, verbose=false)
            bouss_opts = BoussinesqOptions()
            
            t = @elapsed try
                solve_boussinesq_2d(ns_opts, bouss_opts)
            catch
                # Might fail due to time/memory constraints
            end
            return t
        end
        
        # Test small problem sizes
        t1 = time_solve(1)
        t2 = time_solve(2)
        
        @test t1 ≥ 0  # Should complete
        @test t2 ≥ 0  # Should complete
        @test t2 ≥ t1  # Larger problem should take longer (usually)
    end
end

end  # Main testset

# Utility functions for tests

"""Helper function to create a simple temperature field for testing."""
function create_test_temperature_field(nx, ny, T_min=290.0, T_max=310.0)
    x = LinRange(-1, 1, nx)
    y = LinRange(-1, 1, ny)
    
    T = zeros(nx, ny)
    for (i, xi) in enumerate(x), (j, yj) in enumerate(y)
        # Simple temperature distribution
        T[i,j] = T_min + (T_max - T_min) * (1 + xi + yj) / 4
    end
    
    return T, x, y
end

"""Helper function to validate basic result structure."""
function validate_multiphysics_result(result::MultiPhysicsResult)
    # Check required fields exist
    @test !isnothing(result.u)
    @test !isnothing(result.v)  
    @test !isnothing(result.T)
    @test !isnothing(result.x)
    @test !isnothing(result.y)
    
    # Check consistent dimensions
    @test size(result.u) == size(result.v) == size(result.T)
    @test length(result.x) == size(result.u, 1)
    @test length(result.y) == size(result.u, 2)
    
    # Check physical reasonableness
    @test all(isfinite.(result.u))
    @test all(isfinite.(result.v))
    @test all(isfinite.(result.T))
    @test result.T[1,1] > 0  # Temperature should be positive (Kelvin)
    
    # Check result metadata
    @test result.solve_time ≥ 0
    @test result.iterations ≥ 0
    @test !isnothing(result.physics_diagnostics)
    
    return true
end

"""Helper to check if a temperature field has reasonable gradients."""  
function validate_temperature_gradients(T, threshold=1000.0)
    nx, ny = size(T)
    
    # Check that temperature gradients aren't too extreme
    for i in 2:(nx-1), j in 2:(ny-1)
        grad_x = abs(T[i+1,j] - T[i-1,j])
        grad_y = abs(T[i,j+1] - T[i,j-1])
        
        @test grad_x < threshold  # Reasonable gradient
        @test grad_y < threshold
    end
    
    return true
end